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Abstract 

We perform molecular dynamics simulation of water using the TIP5P model to quantify struc- 
tural order in both the first shell (defined by four nearest neighbors) and second shell (defined by 
twelve next-nearest neighbors) of a central water molecule. We find the anomalous decrease of 
orientational order upon compression occurs in both shells, but the anomalous decrease of trans- 
lational order upon compression occurs mainly in the second shell. The decreases of translational 
and orientational orders upon compression ("structural anomaly") are thus correlated only in the 
second shell. Our findings quantitatively confirm the qualitative idea that the thermodynamic, 
dynamic and structural anomalies of water are related to changes in the second shell upon com- 
pression. 
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It is known both from experiments |2j, y, |J] and simulations that the first shell 

of a central water molecule, usually defined by the first minimum of the oxygen-oxygen (O- 
O) pair correlation function (PCF), can accomodate between four and five water molecules, 
depending on pressure [8|. The signature of this first shell, defined by the first maximum of 
the 0-0 PCF, barely changes with pressure. In contrast, the properties of the second shell, 
which extends between the first and second minima of the 0-0 PCF, are highly dependent 
on pressure, indicating that large strnctnra. changes occur in this she.l upon compression fl. 

The structural order of water has been quantified by two measures [9fl: a local orien- 
tational order parameter q, which quantifies the extent to which a molecule and its four 
nearest neighbors adopt a tetrahedral local structure in the first shell, and a translational 
order parameter t, which quantifies the tendency of molecular pairs to adopt preferential 
separations. While q depends only on the four nearest neighbors of a central molecule in its 
first shell, t depends on all the neighbors of the central molecule. 

Water in the liquid phase displays: (i) a thermodynamic anomaly (density decrease upon 
cooling or, equivalently, entropy increase with pressure); (ii) a dynamic anomaly (increase 
of diffusivity upon compression); (iii) a structural anomaly (decrease of both q and t upon 
compression) [9j. Several other liquids with local tetrahedral order 1CJ, III] such as silica, 
silicon, carbon and phosphorous also show waterlike anomalies. In the case of water |9j and 



silica 



12], computer simulation studies show that the anomalies (i)-(iii) in the liquid phase 



are closely related. For example, in the case of water, the region of thermodynamic anomaly 
in the temperature-density (T-p) plane is enclosed by the region of dynamic anomaly, which 
in turn is enclosed by the region of structural anomaly 9j. 

Recent studies show that simple liquids interacting via spherically-symmetric potentials 



can exhibit waterlike anomalies 
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Q, 3' suggesting that strong orientational interac- 
tions in the first shell are not necessary for a liquid to show thermodynamic, dynamic and 
structural anomalies. In light of these findings, it remains unclear how much the strongly 
orientation-dependent first-shell interactions and the second-shell environment each con- 
tribute to water's anomalies. To address these questions, we first modify the definition of 
first and second shells for the purpose of quantitative study. Then we define the orientational 
and translational order parameters within each shell and study their changes with T and p. 

We perform constant volume isothermal (NVT) molecular dynamics simulation of 512 
TIP5P (five-site transferable interaction potential) water molecules. Our simulations are 
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performed using a cubic box with periodic boundary conditions. We control the tempera- 
ture using a Berendsen thermostat [if]]. The TIP5P model reproduces the thermodynamic 
properties of water over a broad region of the phase diagram In particular, we find 

that the TIP5P model reveals similar relations between the thermodynamic, dynamic and 
structural anomalies as observed in ref. [9( using the SPC/E model (see Fig. [T]). 

The first and second shells of water can be defined according to the first and second 
minima of the PCF, g(r). For this definition, the number of molecules in each shell will 
change with pressure and temperature 0, Qj. But the orientational measures that most 
concern us are the tetrahedral arrangement of nearest neighbors, and bond orientational 
order in next-nearest neighbors of a central molecule. To see how these orders evolve across 
a broad range of state points, we must base our comparison on a fixed number of nearest 
and next-nearest neighbors. Moreover, the minima in g(r) become not obvious at high p, 
and g(r) becomes almost featureless beyond the first peak at high p (see Fig. E^b)). Hence 
we choose a less ambiguous shell definition by denoting the nearest four and next-nearest 
twelve neighbors of a central water molecule as the first and second shells respectively. 

We first study the average effect of density on different shells by dividing g(r) into three 
regions. We compute the average number of neighbors of a central molecule at a distance r 
as N(r) = Ann r' 2 g{r')dr' , where n is the number density. We define 7*1 and such that 
iV(ri) = 4 and N(r 2 ) = 16. Therefore, we can define three regions: < r < ri (first shell), 
r\ < r < r2 (second shell), and r > r2, where r\ and r2 depend on T and p. Figure EJ^a) 
shows N(r) at T = 280 K and p = 1.00 g/cm 3 (n = 33.4/nm 3 ), where r\ = 0.32 nm and 
r2 = 0.48 nm. 

Fig. [2](b) shows the 0-0 PCF of TIP5P water at T = 280 K and a range of density 
covering the anomalous regions of water of Fig. [IJ Figure [2](c) shows the change upon 
compression, Ag{r) = g{r)\ p — g(r)\ po , where p = 0.88 g/cm 3 . Figure [2](d) shows the 
corresponding change, AiV(r) = N(r)\ p — N(r)\ Po . Fig. [2](b) shows that as p increases, the 
first peak of g(r) decreases, so Ag(r) < at r = 0.28 nm in Fig. E](c). This effect of p on 
g(r) is a result of having a fixed number of neighbors at r « 0.28 nm, normalized by n in the 
definition of g(r). The change of the number of neighbors corresponding to the first peak of 
g(r) is barely distinguishable (see Fig. E](d)), i.e. AN(r) m — 0.2 for r pa 0.28 nm. This 
implies that the distance defined by the first peak of g(r) is practically impenetrable and 
thus, resembles a hard core. The main changes in g(r) (Fig. EJ^b)) and Ag(r) (Fig. E](c)) 
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occur in the second shell. As the density increases, hydrogen bond bending allows water 
molecules in the second shell to shift toward the first shell, filling the interstitial space ^]. 
The changes of g(r) with pressure for r > r 2 are minimal. 

Figure^d) shows in double logarithmic scale the relationship between AN(r) and r. The 
slope of curve, if), characterizes the power law dependence AN(r) oc r^. There are three main 
regimes in the behavior of AiV(r) as shown by the different slopes if) > 3, if) < 3 and if) = 3. 
The if) = 3 at r > r 2 is mainly caused by the density change, since g(r) ~ 1 (Ag(r) ~ 0) for 
r > r 2 , so AiV(r) behaves approximately as AN(r) oc Ap r 3 with Ap = p — p . Both the 
if) > 3 and if) < 3 regimes are located within the second shell. The increase of AN(r) for 
r < r 2 is not only due to density increase, but also due to the shift of water molecules from 
the second shell around 0.45 nm toward the first shell around 0.28 nm. Thus, the regime 
where if> > 3, for roughly r < 0.4 nm, is due to an increase of g(r) (Ag(r) > 0), and if> < 3 for 
roughly 0.4 nm < r < r 2 is due to the decrease of g(r) (Ag(r) < 0). Note that AN(r) = 1 
for r ~ 0.33 nm, corresponding to the fifth neighbor [8J, which is very close to the border of 
the first shell, where Ag(r) has its maximum value. This fifth neighbor in the vicinity of the 
first shell of water can produce a defect in the tetrahedral network of water at high density. 
This defect leads to hydrogen bond bifurcation and offers paths with low energy barriers 
between different network configurations of water. It is also related to diffusion anomaly by 
lowering energy barriers for translational and rotational motions of water molecules 

The translational order parameter t is defined in refs. jol, Q, 15. 18] 
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t= / \g(s)-l\ds, (1) 
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where the dimensionless variable s = rn 1 ^ 3 is the radial distance r scaled by the mean 
intermolecular distance n _1//3 , and s c usually corresponds to half of the simulation box 
size, which is large enough to have g(s c ) w 1. We can decompose the translational order 
parameter t into ti,t 2 , and £3 for each shell of water by integrating \g(s) — 1| over the three 
different regions < s < Si, Si < s < s 2 and s > s 2 , where s± = rin 1 ^ 3 and s 2 = r 2 n 1 ^ 3 . 
Obviously 

t = tl +t 2 + t 3 . (2) 
The orientational order is used to quantify the tetrahedrality of the first shell, defined 
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as 



-i 2 



cos 9jik + 



(3) 



* = !-§£ £ 

j=i fc=j+i 

^■jfc is the angle formed between neighbors j and k and the central molecule i. The average 
value q = J2iLi 1i quantifies the orientational order of the system based on the molecules 
in the first shell. For perfect tetrahedral order, q — 1; for an uncorrelated (ideal gas) system, 
q = 0. 

Because the second shell of the hexagonal ice crystal forms an hep lattice, the orientational 
order parameter for the second shell of water can be characterized by which quantifies 
the extent to which a molecule % and twelve of its neighbors adopt the local fee, bec, or 
lep structures. This orientational order parameter 19j is often used for simple liquids 
15, Q, 0] because they have fee or bec crystal structures. In order to compute we first 



define twelve bonds connecting each water molecule i with its twelve next-nearest neighbors 
in the second shell, and compute for each bond its azimuthal and polar angles (9, if). Next 
we compute Y em (9, ip), the average of the spherical harmonic function over the 12 bonds of 
the molecule i. Finally we compute 

i 

m=i 



Qt 
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2£ + 



(4) 



For £ = 6, the average value Q 6 = J2i=i Qei quantifies the orientational order of the 
system based on the molecules in the second shell. Qq is large 19| for most crystals such as 
fee (0.574), bec (0.511), hep (0.485). For uncorrelated systems, Q 6 = 1/Vl2 = 0.289. 

Figure [3] shows the density dependence of all six order parameters at three temperatures 
covering the anomalous region of TIP5P water (see Fig. [1]). Although t\ is much larger 
than t2 and £3, it is apparent that £2 makes the most important contribution to the anomaly 
of t (decrease of t with increasing density), compared to t\ and t s . t\ also makes a small 
contribution to the t anomaly at low T = 240 K due to a small decrease in the first peak 
of g(r) upon compression. The anomalous behavior of t becomes weak at T = 280 K and 
disappears at T = 320 K. The orientational order parameters q and Qq both show similar 
anomalous behavior. The distribution of individual qi shifts from high q (ice-like) at low 
p and T to low q (less tetrahedral) at high p and T as shown in Fig. H](a) and (d), due to 
increased hydrogen bond bifurcation [8] as interstitial molecules move closely to the first 
shell (Fig. Hfc) and (f)). always has approximately normal distribution as shown in 



Fig. H](b) and (e) because there is no direct bonding between center water molecule and 
second shell of water. 

A useful way of investigating structural order in fluids is to map state points onto the t-q 
plane, a representation called the order map 0, [20]. The order map for TIP5P water 
(i.e., using t and q) is shown in Fig. E^a). This order map is similar to the one found in 
ref. js| using the SPC/E model. Its main characteristic is the correlation of the two order 
parameters in the anomalous regions where both q and t decrease with density, as shown by 
the isotherms collapsing onto a line. Fig. H^b)-(h) shows the different order maps obtained 
by considering the order parameters in different shells. The only one that shows the states 
in the thermodynamically, dynamically and structurally anomalous regions collapsing onto 
a line, is the panel (f) (i.e. the t 2 -Qe order map of the second shell), indicating that the 
changes in the second shell are related to anomalies of water. 

The first shell order map t\-q in (c) is not correlated because t\ has only small changes 
with increasing density due to the impenetrable hard core at 0.28 nm, while q changes 
significantly with density. In the second shell, £2 and Qq both change significantly and 
simultaneously with density so that they are well correlated. Our work quantitatively shows 
that the second shell is related to anomalies of water by its gradual shift towards first shell 
upon compression. In addition to water, other tetrahedral liquids such as silica, silicon, 



carbon and phosphorous 



lOj may also exhibit similar behavior, and a detailed, shell-based 



study of their order parameters may prove useful. 
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FIG. 1: (Color online) Three anomalous regions in T-p plane for 512 molecules interacting with 
TIP5P potential, (i) The density anomaly region is defined by the locus of density maxima (TMD), 
inside of which the density increases when the system is heated at constant pressure, (ii) The 
diffusion anomaly region is defined by the loci of diffusion maxima or minima (DM), inside which 
the diffusivity increases with density, (iii) The structural anomaly region is defined by the loci of 
translational order minima (i m i n ) and maxima (t m ax)j or orientational order maxima (9max)> inside 
which both translational and orientational orders decrease with density (see Fig. [3]) . 
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FIG. 2: (a) The number of neighbors N(r) around a central water molecule. r\ and r<i correspond 
to the first and the second shell distances, defined such that N(r±) = 4 and N(r2) = 16. (b) 
The 0-0 PCF g(r), (c) Difference Ag(r) between g(r) at a given density and g(r) at po, and (d) 
Difference AiV(r) between N(r) at a given density and N(r) at po for TIP5P water, ip characterizes 
the local slope. The bold portions of the curves correspond to water's second shell, r\ < r < r2, 
showing that the largest changes upon compression occur in the second shell. 
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FIG. 3: Translational order parameters t (total), t\ (first shell), ti (second shell), ts and orienta- 
tional order parameters q (first shell), Qq (second shell) of TIP5P water as function of density at 
different T. The anomalous decrease of orientational order upon compression occurs in both shells 
{QiQs)i but the anomalous decrease of translational order upon compression mainly occurs in the 
second shell (<2). 
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FIG. 4: Histograms of (a) the local orientational order qi in the first shell, (b) in the second 
shell, and (c) distance r§i between a central water molecule i and its fifth neighbor of TIP5P water, 
(a), (b), and (c) show the changes for three different p at fixed T = 280K. (d), (e), and (f) show 
the changes for three different T at fixed p = 1.00 g/cm 3 . Upon compression or heating over 
anomalous regions of phase diagram, the fifth neighbor (and other interstitial water molecules in 
the second neighbor shell) shift towards first shell (see also Fig. [2] and ref. corresponding to 
anomalous changes of structural order in the first and second shells as quantified by Fig. [3] and 
Fig. [3 



11 



1.5 



0.8 











r(a)' ' ' ' y 

:\\^ 


I 1 1 1 1 1 1 >■ 

( b ) : 

L/jT ■ T=240KJ 

• ■ T=280K- 

T=320KT 


r(c) 


(d) -i 






r(e) 


(f) -. 


: (g) 

s- * ' - - ■ ■ ' " 


(h) : 


I . I . I 





0.5 0.6 0.7 0.8 0.26 0.27 0.28 0.29 0.3 

q Q« 



FIG. 5: Order maps for TIP5P water (color online). The arrows indicate the direction of increasing 
density from 0.84 g/cm 3 to 1.32 g/cm 3 . Only for the second shell order map, t2-Q% in (f), the 
isotherms collapse on a line and the decrease of translational and orientational orders is correlated. 
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